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We study chaotic behavior of order parameters in two coupled ensembles of self- 
sustained oscillators. Coupling within each of these ensembles is switched on and off 
alternately, while the mutual interaction between these two subsystems is arranged 
through quadratic nonlinear coupling. We show numerically that in the course of 
alternating Kuramoto transitions to synchrony and back to asynchrony, the exchange 
of excitations between two subpopulations proceeds in such a way that their collective 
phases are governed by an expanding circle map similar to the Bernoulli map. We 
perform the Lyapunov analysis of the dynamics and discuss finite-size effects. 
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Behavior of high-dimensional nonhnear systems of different nature (e.g. in hy- 
drodynamics, electronics, neurodynamics, laser physics and nonlinear optics) in 
many cases can be treated in terms of a cooperative action of a large number of 
relatively simple elements, such as interacting oscillators. Quite often the essen- 
tial dynamics is low-dimensional: even for a very large number of elements (and 
also in the thermodynamic limit where this number tends to infinity) one can 
find a few degrees of freedom that describe the macroscopic evolution. A famous 
example is the Kuramoto model where the global variable — complex order pa- 
rameter — undergoes a Hopf bifurcation corresponding to self-synchronization in 
the ensemble of oscillators. Here we report on a modification of the Kuramoto 
model where the behavior of the global variables becomes chaotic. Our model 
consists of two populations that undergo Kuramoto-type transitions, where due 
to an external modulation of the coupling strength the oscillators synchronize 
and desynchronize alternately in both subpopulations. The operation of the 
whole system consists of alternating patches of activity of the subensembles, in 
the sense that their order parameters (complex mean fields) vary between small 
and large magnitudes. Due to a specially constructed additional nonlinear cou- 
pling between the subpopulations, the dynamics of the phases of the complex 
order parameters is described stroboscopically by an expanding circle map (the 
Bernoulli map). Apparently, such systems may be designed on a base of ensem- 
bles of electronic devices, like arrays of Josephson junctions, or with nonlinear 
optical systems, such as arrays of semiconductor lasers. 



I. INTRODUCTION 

A long-time challenging problem motivating development of nonlinear science concerns 
complex behavior of systems characterized by a large number of degrees of freedom. Such 
systems are of interest in various fields, e.g. in hydrodynamics (the problem of turbulence), 
in laser physics and nonlinear optics, electronics, neurodynamics^. The functioning of such 
systems often can be thought of as cooperative action of a large number of relatively simple 
elements, such as oscillators, with different types of interaction between them^i"-. In many 
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cases, nevertheless, the dynamics can be effectively treated as low- dimensional in terms of 
appropriate collective modes^^"-'^"-. In this context, an interesting question arises about 
a possibility of implementing various known types of low- dimensional complex dynamic 
behavior on the level of the collective modes. 

One basic textbook example of strong chaos is Bernoulli map, or expanding circle map^. 
In general form, the map is v'n+i = Mipn (mod 27r), where M is an integer larger than 1, 
and is an angular variable. With initial condition for ip/2Ti represented in the numeral 
system of base M by a random sequence of digits, the dynamics will correspond to a shift 
of this sequence by one position to the left on each next step of the iterations. It means 
that the representative point will visit in random manner M equal segments partitioning 
the circle, in accordance with the mentioned sequence of the digits. This is just what we 
call chaos. The sensitive dependence of the dynamics on initial conditions is characterized 
quantitatively by a positive Lyapunov exponent A = In M. 

In a recent series of papers^""— it was shown how to implement the dynamics described by 
the Bernoulli map in low-dimensional systems of alternately excited self- sustained oscillators. 
In these systems the existence of uniformly hyperbolic chaotic attractors of Smale - Williams 
type (see, e.g., Ref.— ) was establishedr^'^^. 

Chaotic nature of the dynamics reveals itself in the chaotic evolution of the phases of 
oscillations generated at successive stages of excitation of the subsystems. The purpose of 
the present article is to demonstrate a possibility of chaotic behavior of similar nature at 
the level of collective variables for multidimensional systems represented by ensembles of 
oscillators. 

We do not start here with some arrangement motivated by a particular application, but 
construct an idealized example that is simple and convenient for the theoretical description. 
The model is composed of two similar ensembles of self-sustained oscillators with their 
natural frequencies distributed in some range; within each of these ensembles the global 
coupling is switched on and off alternately. Interaction between these two subensembles is 
arranged with a special type of additional coupling through mean fields characterized by 
quadratic nonlinearity. 

Due to presence of the global coupling, each ensemble can undergo the Kuramoto tran- 
sition: as oscillators synchronize, the collective field emerges with notable amplitude and 
definite phase of its oscillations. As the coupling is turned off, the oscillators desynchronize 
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because of the frequency detuning of individual oscillators, and the collective field disap- 
pears. The operation of the whole system consists of alternating activities of two ensembles 
with the corresponding alternating meandering of their order parameters. Due to an addi- 
tional coupling, the excitation is transmitted from one ensemble to another, and back, so 
that the phase of this collective excitation in the course of the process evolves in accordance 
with the expanding circle map mentioned above. 

One of the goals of this paper is to compare the global realization of hyperbolic chaos, 
as outlined above, with other types of collective chaos in ensembles of dynamical systems. 
Postponing a detailed discussion to the conclusion, we mention here that collective chaos 
has been mainly studied in two different contexts: in ensembles of maps^"— , where indi- 
vidual elements are chaotic, and in populations of oscillators (which as individual elements 
are, contrary to the case of maps, nonchaotic), either with a distribution of natural frequen- 
pjggi9i2o ifigntical^'Si"— . In this context our study is closely related to that in^^, because 
we also consider non-identical oscillators; the difference is that we organize the coupling in 
a special way to ensure desired properties of collective chaos. 

II. BASIC MODELS 

In this Section we formulate models of interacting oscillator ensembles, using three dif- 
ferent levels of reduction. First, we consider the equations in an "original" form, where each 
oscillator is described by the van der Pol equation. Next, we exploit the method of averag- 
ing and formulate the model in terms of slowly varying complex amplitudes of oscillations. 
Finally, we proceed with neglecting amplitude variations for single oscillators, and account 
only for variations of the phases on their limit cycles; that is the model of ensemble of phase 
oscillators. 

As outlined in the Introduction, we do not study a general system of coupled alternately 
synchronized ensembles, but construct a model that should produce chaos, having specific 
properties (hyperbolicity). Moreover, we want to check, how reductions to amplitude and 
phase equations influence the dynamical properties. Therefore the model in "original" vari- 
ables (Eq. ([3])) below) looks rather cumbersome, while its reduction to complex amplitude 
(Eq. ([5])) and phase (Eq. ([8])) variables are rather close to the standard Kuramoto model. 
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A. Coupled van der Pol oscillators 

Let us consider two interacting ensembles of self-sustained oscillators. We assume that 
the sizes K of these ensembles are the same, and the oscillators are described by variables 
and Uk, respectively, where k — 1, . . . , K . Next, we assume that the distributions of natural 
frequencies cuk of the oscillators are identical for both ensembles, with some mean frequency 
cuq. Within each ensemble, oscillators are coupled via their mean fields X and Y, defined 
according to 

k k 

To account for the dissipative nature of the coupling (i.e. a tendency to equalize the in- 
stant states of the interacting subsystems), we assume that it is introduced by terms in the 
equations containing time derivatives of the fields X and Y. We suppose that the coupling 
strength varies in time periodically, with a slow period T ^ Stt/o'o, between zero and some 
maximum value, which exceeds the synchronization threshold of the Kuramoto transition. 
Thus, each ensemble goes periodically through the stages of synchrony (when the coupling 
is large) and asynchrony (when the coupling is small). The coupling is organized in such 
way that these stages in the two ensembles occur alternately. Finally, there is a nonlinear 
interaction between ensembles via the second-order mean fields 

k k 

Being represented as sums of the squares of the original variables, these fields X2 and 
Y2 contain components with the double frequency in comparison to the frequency of the 
variables X and Y. To ensure an efficient resonant interaction we need the components 
with the main oscillator frequency; to this end the coupling terms are chosen as products of 
time derivatives of the fields ^2,12 and of an auxiliary signal sinouot. These products then 
contain the components with the basic frequency ujq. The set of governing equations for the 
model reads: 

dt^ 

d^Vk 
dt^ 

where A; = 1, 2, X. The functions fi{t) = cos^ (nt/T) and f2{t) = sin^(7ri/T) describe the 
alternate on/off switching of the couplings inside the ensembles. Parameter Q determines the 
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(3) 



amplitude of each single van der Pol oscillator; parameters k and e characterize, respectively, 
the internal and mutual couplings for the ensembles. It is assumed that the period of 
coupling modulation contains a large integer number of periods of the auxiliary signal, i.e. 
WoT/27r = » 1. 



B. Description in terms of complex amplitudes 

The method of averaging allows us to describe weakly nonlinear oscillators in a simplified 
form, via the slow dynamics of the complex amplitudes; this description is appropriate under 
assumptions that \ujk — ujq\ « wq, 2ti /T « uq. For the van der Pol oscillators ([3]) the 
averaging can be accomplished by introducing complex amplitudes a^, bk according to 

a,{t) = e-^^otit±^ , b,{t) = e-^^y^±^ . (4) 
iujo iojQ 

Substituting this in Eqs. ([3]) and averaging over the period of fast oscillations 271 /ojq (prac- 
tically, this may be done simply by dropping all term in the r.h.s., which contain fast time 
dependencies like ~ e^*'^"*, e'*^^*'^"* etc.), we obtain 

dfc = iVlkCik + ^(<5 - la/cHcifc + ^Kfi{t)A + ^^£^^2 , 

(5) 

hk = iVtkbk + -(Q - l&fcp) bk + -Kf2{t)B + -ieA2 . 

Here Qk = — are the frequency differences with respect to the average one. The mean 
fields A, B, A2, B2 are defined similarly to ([T]) and ([2]) as 

•4=^E''». ^^^E"- -4, = 15:4, B. = ^5:(.|. (6) 

k k k k 

C. Phase approximation 

One more step in simplification of the model is to completely neglect the amplitude 
variations for single elements, and to reduce the description to that in terms of ensembles 
of phase oscillators. In this approximation the individual oscillators are assumed to have 
a constant amplitude (that of the limit cycle of a single oscillator) throughout the process, 
while the dynamics manifests itself only in the evolution of their phases. This is a widely 
used approximation in the theory of synchronization, see^i^. 



To derive the phase equations, we substitute in Eq. ([5]) 

afc(t) ^ i?e'^'=W , bk{t) ^ Re'^"^'^ , (7) 

where is a constant, equal for all oscillators. We specify R = \/Q + k,/2 according to the 
following reason. First, if one switches off the coupling between the ensembles setting e = 0, 
then each single oscillator is described by the equation a = iVLa + ^{Q — \a\^)a + \k,A. In 
absence of the synchronization A = 0, and the amplitude of the stationary self-sustained 
oscillations is determined from |ap = Q. On the other hand, in the case of maximal 
synchronization one has / = 0, A = a, and, hence, |ap = Q + Since the dynamics 
consists of the alternating epochs of synchronization and desynchronization, it is reasonable 
to take the average value.— 

Then, from Eqs. ([5]), we obtain equations for the phase dynamics: 

eu = nu + hm[{t,h{t)U + ieRV2)e-''^] , 

\ (8) 
0fc = fi, + -Im[(/€/2(t)V + z£i?t/2)e-^<^^] , 

where the complex mean fields [/, V, U21 V2 are defined according to 

k k k k 

D. Thermodynamic limit 

Above we have assumed the number of oscillators in the ensembles to be finite. Now let 
us write the equations in the "thermodynamic limit" K ^ 00. In the case we deal with, the 
oscillators differ only by their natural frequencies, then it is convenient to parameterize the 
oscillators by this continuous variable, i.e. by the frequency cj, and introduce a distribution 
over the frequencies characterized by a function g{uj). Respectively, after transformation of 
the equations to slow amplitudes, and in the phase approximation, we will use the index 
VL = oj — ojq designating the frequency difference. Then, we characterize the distribution by 
the density giVt) = g{ujQ + Vt). 

With these notations, the set of the phase equations (IE]) can be rewritten as 

9n = n + ^Im [{Kh{t)U + ieRV2)e~''''] , 

1 (10) 
0n = fi + -Im[(/€/2(t)V + ^£i?f/2)e-^'^"] , 
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where the mean fields are defined now via the integrals over the density 



r r (11) 

U2= dn g{n)e^'^» , V2= dn ~g{n)e^"^^ . 

In a similar way, the equations for the complex amplitudes and for the original van der Pol 
oscillators may be easily reformulated in the thermodynamic limit. 



III. NUMERICAL EVIDENCE FOR COLLECTIVE CHAOS 

In this Section we present numerical studies of the models (jS]), ([5]), and ([8]). 

As is known from the theory of synchronization in populations of oscillators developed by 
Kuramoto, the properties of the synchronization transition are qualitatively the same for all 
unimodal smooth distributions of oscillators over their frequencies. In computations below 
we specify the distribution for the model ([3]) as 

(TT -niuO — UJq — Aw) r a at 

-COS ^ , a;G[a;o-Ac.,c.o + Aa;], ^^^^ 
otherwise , 

and set ujq = 27r and Au = |. 

We select the form (fT2l) because it allows us to choose the discrete set of frequencies for 
a finite ensemble according to a simple relation 

'2k - 1 



2Au 

Uk = ujq — Aw H arccos 

TT 



- 1 



1,2,. ..K. (13) 



K 

Furthermore, absence of significant tails of the distribution (compared, e.g., to a Lorentzian 
one) simplifies numerical studies, as we do not have to bother about non-resonant oscillators 
with too large or too small frequencies. For the models ([5]), and ([8]) analogous distributions 
were used obtained from f|T2l) and f|T3l) with the substitution Q = u — Uq. 

With this setup, we simulated dynamics of the models (jS]), ([5]), and ([8]) in computations at 
parameters T = 100, k = 1, e = 0.1, and Q = 3, for system sizes K = 1000 and K = 10000. 
The main quantities of interest are the phases of the mean fields. For the ensembles of van 
der Pol oscillators we define them according to 

X Y 

$x = ~ arctan , $y = — arctan (14) 

OJqX OJqY 
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For the ensembles ([5]) and ([8]) we define the phases simply as the arguments of the complex 
mean fields 



A = lAle^*^ , B = \B\e' 



U=\U\e'^^, V=\V\e 



(15) 
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FIG. 1. The mean fields (panels (a) and (b)) and their phases (panel (c)) for two interacting 
ensembles of van der Pol oscillators (jl2p with K = 1000. The modulation function of the coupling 
/i^2 are shown in (a,b) with dashed lines. To make the picture (c) clear, we subtract a linear phase 
growth and plot the value ^x,Y — wi, where u = ujq — 0.2 is chosen close to the empirical mean 
frequency u , which slightly differs from loq due to the nonlinearity. — is shown with solid 
red line and $y — ut with dashed blue one. 



Figure [T] illustrates the dynamics of the mean fields X{t) and Y{t) in the ensemble of 
van der Pol oscillators ([3]). First, we mention that because of the modulation ~ f^fi,2{t) the 
mean fields of two ensembles vary significantly: they drop nearly to zero in the epochs where 
the corresponding values of /i,2 are small, and attain large values in the epochs where the 
coupling inside a population becomes large. In panel (c) we show the time dependences of the 
phases of the mean fields $x.y • One can see slight regular variations of the phases correlated 
to the amplitudes of the mean fields, and additional phase shifts close to the moments of 



time when the corresponding mean fields nearly vanish (at times t ~ 40, 140, 240 for X{t), 
and at times t ^ 90,190,290 for Y{t)). These phase shifts correspond to transfer of the 
phase from one ensemble to another through their mutual coupling terms proportional to e. 

Usually, as an ensemble of oscillators passes through the Kuramoto transition from a 
non-synchronous to a synchronous state while the coupling strength increases, the phase 
(potentially, of arbitrary value) of the arising collective mode is determined by fiuctuations 
stimulating the excitation of this mode. In our setup, however, the excitation occurs in 
presence of a small driving force ~ e because of action of another ensemble, which is syn- 
chronous at that moment and generates notable mean field. This stimulation determines 
the phase on the appearing collective mode, which accepts this externally designated phase. 
This is the mechanism of the phase transfer. 

Because the mutual couplings are proportional to the second-order mean fields character- 
ized by a doubled frequency, the phase transfer is accompanied with doubling of the phase. 
To explain this, let us assume that at the transfer of excitation from the first to the second 
subensemble we have X ~ cos(a;ot + $) and, respectively, X2 ~ sin[2(ci;ot + $) + const]. 
Then, the driving force affecting the second subensemble contains the resonance component 
of X2 sin Wot ~ cos(u;o'^ + 2$ + const). Respectively, the arising mean field will be of the 
form Y ~ cos(a;ot + 2$ + const). 

As doubling of the phase occurs at each transfer from one subsystem to another and 
back through the full cycle (i.e. over the period T), as a result we expect that the phase is 
multiplied by the factor of 4 (up to an additive constant). 

To check the supposed mechanism of the phase transfer, we constructed numerically a 
stroboscopic iteration phase diagram ^x{nT) $x((^^ + 1)^), relating the phases at the 
successive moments of maximal amplitude of the mean field X. This map, shown in Fig.[3]^a), 
clearly indicates that the transformation is indeed close to the Bernoulli-type map 

^x{{n + 1)T) = A^x{nT) + const (mod 27r) . (16) 

For the same values of the parameters we also simulated the dynamics of the coupled 
oscillator ensembles in the slow complex amphtude version of Eq. (E]) (see Fig. |21^a,b)), and 
in the phase approximation of Eq. ([8]) (see Fig. |^c,d)). In these situations the second-order 
effects of amplitude-dependent frequency shifts, like those seen in Fig. are not observed. 
As a result, the phases between the short intervals of phase transfers are nearly constant. 
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and the phase shifts at the transfers are clearly visible. The stroboscopic maps of the phases 
are shown in Fig. |3t^b,c). They demonstrate Bernoulli-type maps similar to that for the 
ensemble of van der Pol oscillators of Eq. (ITBil (see Fig. [3](a)). 




200 400 600 800 1000 200 400 600 800 1000 



time t time t 



FIG. 2. Evolution of amplitudes and phases for the mean fields of coupled ensembles described by 
slow complex amplitudes (panels (a,b)) and in the phase approximation (panels (c,d)), K = 1000. 
Solid lines: variables |A|,<I>A; \U\,^u, dashed lines: variables \B\,^b, \V\,^v- 

For a further numerical characterization of the chaotic dynamics of mean fields, we con- 
structed the stroboscopic maps of the complex mean fields via period of modulation T. 
Portraits of attractors of these maps for our three levels of description of the ensembles are 
depicted in Fig. HI For a dynamical system where the phase (or another cyclic variable) 
undergoes a Bernoulli-type transformation, while in other directions in the phase space the 
phase volume compresses, one expects the strange attractor to be of the Smale- Williams 
type, i.e. represented by a solenoid. In a two-dimensional projection this attractor looks 
like a circle with a fractal transversal structure. Fig. H] confirms this picture for the dynamics 
of the mean field. 

We would like to stress that the Bernoulli map describes the collective phase (i.e. the 
phase of the mean field) but not individual phases of the oscillators. Indeed, as is evident 
from topological considerations, the doubling of the phase is only possible, if the amplitude 
vanishes at some stage. For the complex mean field the phase doubling is achieved by 
synchronization - desynchronization, while all individual oscillators always have a finite 
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FIG. 3. Stroboscopic maps over the period of external modulation T for: (a) the ensembles of 
van der Pol oscillators p^ : (b) for the ensembles described by slow complex amplitudes ([5]); 
(c) for the ensembles in the phase approximation ([8]); in all cases K = 1000. In all cases the 
dynamics seems well described by the Bernoulli map (|16p . The observed splitting of the "lines" 
(the most pronounced in panel (a)) appears because of presence of transversal fractal structure of 
the attractors (see FigU]): distinct filaments of the attractor give rise do distinct filaments on the 
phase iteration diagram due to imperfection of the phase definition. 




2 -0.7 0.7 




0.58 
0.58 



FIG. 4. Projections of the stroboscopic maps on the plane of the order parameters: {X,X) for 
the van der Pol oscillators (left column); (Reyl, ImA) for the ensembles described by slow complex 
amplitudes (center column); (ReC/, Im{7) for the ensembles in the phase approximation. Bottom 
row show enlargements to make the fractal transversal structure evident. Here K = 10000. 
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amplitude and therefore cannot be described by the doubhng Bernoulh map. We illustrate 
this in Fig. [5l where the maps are shown similar to that of Fig. |3l but for the individual 
phases of two representative oscillators of the ensemble. Mostly close to the Bernoulli map 
is the behavior of the oscillator at the center of the frequency band, as this oscillator nearly 
perfectly follows the mean field. Nevertheless, one can clearly see the "defects" of the 
transformation appearing because the "amplitude" of the oscillator cannot vanish. The 
oscillator at the edge of the band does not generally follow the phase of the mean field, and 
its dynamics is far from the Bernoulli map. 




vr 27r tt 27r 



FIG. 5. Stroboscopic maps over the period of external modulation T for individual oscillators. Left 
panel: oscillator in the center of the band with uj = ujq; right panel: oscillator with u = ujq — Au. 

IV. FINITE SIZE EFFECTS 

In this section we characterize collective chaos for different values of the ensemble size K. 
We focus here on the properties of the model Eq. (jS]) in terms of complex amplitudes. First, 
we analyzed the stroboscopically observed phases of the mean field and found that for some 
ensemble sizes K a periodic behavior is observed. The periods of the found periodic regimes 
are shown in Fig. M^a) vs K. The values of K for which no period is plotted correspond to 
the chaotic states. At first glance, this contradicts robustness of chaos expected because of 
its approximate description in terms of the expanding Bernoulli map. Apparently, changing 
the ensemble size we make an essential perturbation of the effective collective dynamics, so 
that the usual arguments of structural stability do not apply. 

The analysis via the bifurcation diagrams Fig. M^a) is confirmed by calculation of the 
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FIG. 6. Bifurcation diagrams for the ensembles in the complex amplitude formulation ([5]) showing 
dynamical regimes in dependence on the ensemble size K. (a): Periods for periodic regimes. 
The values of K for which no period is plotted correspond to the chaotic states, (b): Three 
largest Lyapunov exponents (the first: black filled circles, the second: red pluses, the third: green 
diamonds). Only a few regimes with small number of oscillators have two positive exponents. In 
panel (c) we show only positive largest Lyapunov exponents, in a logarithmic scale, to demonstrate 
that it does not tend to decrease for large ensemble sizes K. 

Lyapunov exponents for the ensembles. We have performed this analysis for the ensembles 
described by complex amplitudes ([5]) of different size. First, in Fig. [7] we present the full 
spectrum of the Lyapunov exponents for three ensemble sizes. For these parameters we 
observe chaos, and in all cases only one Lyapunov exponent is positive. We interpret this 
as existence of one collective chaotic mode, contrary to the cases when many Lyapunov 
exponents in populations of oscillators are positive (cf.-'^'^). 

In Fig. El^b) we present the three first Lyapunov exponents for the same data as in 
Fig. |6](a). Again, like in Fig. [6](a), one can see that for certain system sizes the dynamics 
is regular, as the largest exponent is negative. One can notice that the largest positive 
Lyapunov exponent decreases with K ioi K < 150, but, as panel Fig. [6](c) shows, saturates 
and presumably does not further decrease for larger system sizes (contrary, e.g. to the 
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FIG. 7. Lyapunov exponents Aj plotted vs index ijK for the ensembles described by complex 
amplitudes ([5]) with K = 20, 30, 100. 2K exponents corresponding to the phases are close to zero, 
while the rest 2K exponents corresponding to the stable amplitudes are negative. Right panel 
shows the region with small in all cases only the first exponent is positive. 

reported dependence \max ~ -ft' ^ in the standard Kuramoto model), although periodic 
windows can be observed for large K as well. (Preliminary calculations demonstrate that 
for large K these windows become extremely rare.) 

We suggest that the observed peculiarities are specific just to finite-size effects, and that 
they will possibly disappear in the thermodynamics limit. Indeed, the computations show 
that they become less expressed under increase of the ensemble size; however, as may be 
concluded from computations, the convergence to the large-size behavior is slow enough. 



V. CONCLUSION 



We have proposed and studied a model system, which demonstrates chaotic dynamics of 
the phases of the mean fields. These collective variables are described by an expanding circle 
map transformation in the course of the transfer of collective excitation alternately between 
two synchronizing and desynchronizing groups of oscillators. We have demonstrated this 
on different levels of description: for the original system of coupled van der Pol oscillators, 
for the model based on the equations for slowly varying complex amplitude, and for the 
phase oscillators. While collective chaos is observed in a wide range of parameters, it is 
not completely robust with the respect to variations of the ensemble sizes: here we observe 
regularity windows. These finite size effects require further investigations. 

An interesting property of the model studied is that in the chaotic state it has only one 
positive Lyapunov exponent (except for several regimes with a very low number of oscillators 
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in ensemble), all others are negative. In this respect our model differs from the ensemble 
of identical oscillators studied in^iSi where the number of positive Lyapunov exponents in 
the regime of collective chaos was macroscopic (proportional to the size of the ensemble). It 
also differs the ensemble of identical Josephson junctions^, where few Lyapunov exponents 
are positive but the macroscopic majority of them vanish due to partial integrability of 
the system. At the moment it is not clear, which physical properties of the systems are 
responsible for this difference. Possibly a comparison with the Lyapunov spectrum in models 
with distribution of frequencies^iSS could shed light on this problem. A promising approach 
for future studies is a calculation of finite-size Lyapunov exponents like in^^. 

Another peculiar property of the system studied is its sensitivity to the number of oscilla- 
tors in the ensemble and appearance of periodic "windows" in dependence on this parameter. 
Such a sensitivity has not been reported for other models demonstrating collective chaos. 
We speculate that this property might be related to the mentioned above existence of one 
positive Lyapunov exponent, what makes the collective chaos in the ensemble less robust. 
This issue requires special attention in the future work. 

We believe, that the model proposed, although rather artificial to be observed in natural 
oscillator ensembles, provides a useful test system for analysis. On the other hand, our 
research opens a possibility of constructing realistic systems with collective chaotic phase 
dynamics based on ensembles of such individual elements that show only regular dynamics. 
This might be feasible, e.g., on the basis of electronic devices, such as arrays of Josephson 
junctions^^, or with nonlinear optical systems, such as arrays of semiconductor lasers^^. Such 
systems are expected to generate robust chaos providing the power level much higher than 
that characteristic for the individual elements. Systems of this kind may be of interest for 
applications requiring generation of chaotic signals, such as communication schemes^S*^, 
noise radar—, etc. 
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